Submitted to Astrophysical Journal Letters 

Preprint typeset using LM^X style cmulatcapj v. 6/22/04 



CONSTRAINED COSMOLOGICAL SIMULATIONS OF DARK MATTER HALOS 

Emilio Romano-Diaz 1,2 , Andreas Faltenbacher 3 , Daniel Jones 2-4 , Clayton Heller 4 , 
Yehuda Hoffman 1 , and Isaac Shlosman 2 

Submitted to Astrophysical Journal Letters 



in 
o 
o 

(N 
bJQ' 

< 



> 
(N 

oo 
o 
in 
o 

^ ' 

Or 

6 ■ 

a: 



ABSTRACT 

The formation and structure of dark matter (DM) halos is studied by means of constrained real- 
izations of Gaussian fields using A-body simulations. A series of experiments of the formation of a 
1O 12 /i _1 M halo is designed to study the dependence of the density profile on its merging history. 
We confirm that the halo growth consists of violent and quiescent phases, with the density well ap- 
proximated by the Navarro-Frenk- White (NFW) profile during the latter phases. We find that (1) the 
NFW scale radius R s stays constant during the quiescent phase and grows abruptly during the violent 
one. In contrast, the virial radius grows linearly during the quiescent and abruptly during the violent 
phases. (2) The central density stays unchanged during the quiescent phase while dropping abruptly 
during the violent phase. (3) The value of R s reflects the violent merging history of the halo, and 
depends on the number of violent events and their fractional magnitudes, independent of the time and 
order of these events. It does not reflect the formation time of the halo. (4) The fractional change in 
R s is a nonlinear function of the fractional absorbed kinetic energy within R s in a violent event. 
Subject headings: cosmology: dark matter — galaxies: evolution — galaxies: formation — galaxies: 
halos — galaxies: interactions — galaxies: kinematics and dynamics 



1. INTRODUCTION 

The problem of the formation and structure of dark 
matter halos constitutes one of the outstanding chal- 
lenges of modern cosmogony and structure formation 
models. The problem is easily formulated as what is the 
outcome of the collapse and virialization of bound per- 
turbations in an otherwise homogenous and isotropic ex- 
panding universe. This is further simplified by consider- 
ing only collisionless non-interacting particles, hereafter 
referred to as dark matter (DM). The resulting collapsed 
objects are dubbed here as halos. This seemingly simple 
problem does not easily yield itself to an analytical un- 
derstanding — the problem addressed here in a series of 
numerical 'experiments'. 

Analytical and numerical studies of the collapse and 
virialization of structures in an expanding universe date 
to the early times of modern cosmogony. The only rigor- 
ous and exact analytical solution relevant to the problem 
is that of a single scale free spherical density perturba- 
tion in a Fri edmann universe, the so-called secondary 
infall model llGunn 1119771 iFillmore fc Goldreich I IT981 
iBertschi nger ||198 5|) and its applic ation to cosmological 
models (Hoffm an fc Shaham 1119851 etc.). The shortcom- 
ing of the analytical approach prompted the study of the 
formation of D M halos by means of N-body simulations 
l|White Ill97fil etc.). The structure of DM halos inferred 
from a variety of cosmological models and power spectra 
of the primordial perturbation field was found to be well 
approximated by a spherically averaged density profile 
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where p s is the density at the scale radius R s . The 
NFW profile constitutes a universal fit to the struc- 
ture of DM halos over many orders of magnitudes in 
mass and at different redshifts. It has been confirmed 
by numerous N-body simulations, although some modifi- 
catio ns have been suggested (most notably iMoore et alJ 
Il998|) . The origin of the NFW profile has been stud- 
ied in the framework of the secondary infall model (e.g., 
Nusser fc ShetjTllT999l ILokas fc Hoffn^l200flt iNTrsserl 



20011 I Asc asil3ajet"ldn I2004J) and the merger scenario 
(jEl-Zantll2005|) . While analytical models necessarily 
invoke spherical symmetry, simulations emphasize the 
background asymmetry. This has led us to embark on 
a series of numerical experiments carefully designed to 
shed light on the problem. 

In this Letter we address the role of the merg- 
ing history in shaping up the density profile. It 
has been determined that the evolution of DM ha- 
los p roceeds in two phases , of rapid and slow accre- 
tion ifWechsler et al. I l200l iZentner fc Bullock 1 12001 
iZhao et al. 1 12003 iSalvador-Sole et alJ 120051 and refs. 
therein). The general understanding that has followed 
is that an NFW structure is quickly established after 
the rapid phase and is preserved during the slow accre- 
tion. Thus it follows that the emergence and evolution 
of the NFW pr ofile might depe nd primarily on the merg- 
ers epoch (e.g.. lEl-Zant ll2005l) . This motivates our first 
study of the formation of DM halos by using constrained 
simulations. We design a set of numerical experiments in 
which a given halo of mass 10 12 /i _1 Mq (where h is the 
Hubble's constant in units of 100 km s _1 Mpc _1 ) is con- 
strained to follow different merging histories. We then 
study the different constrained halos and compare the 
evolution of their density profiles. 
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Constrained Cosmological Simulations of DM Halos 



The ability to design initial conditions for TV-body sim- 
ulations by constrained realizations of Gaussian fields is 
the key to the 'experimen tal' approach used here. T his 
is done by following the iHoffman fc Ribakl 1)1991(1 al- 
gorithm. The design of an experiment starts with the 
expression of whatever constraints one would like to im- 
pose into linear constraints on the primordial perturba- 
tion field, which is assumed to be Gaussian within the in- 
flation paradigm. The constrained realizatio n is used to 
set up the initial conditions by means of the iZel'dovich I 
l|1970() approximation. The resulting simulation is re- 
ferred here as a constrained simulation. 

The linear constraints are described in the models 
are presented in 33 and their analysis in 21 The cosmo- 
logical implications are discussed in 

2. CONSTRAINED SIMULATIONS 

We used the updated version of the FTM-4.4 hy- 
brid code (Heller & Shlosman 1994; Heller 1995), with 
N ~ 2.1 x 10 6 . The gravitational forces are computed us- 
ing the routine f alcDN (Dehncn 2002), which is about ten 
times faster than optimally coded Barnes & Hut (1986) 
tree code. The gravitational softening is e = 500 pc. 
The addressed issue of the collapse of an individual halo 
in an expanding Friedmann universe has led us to use 
vacuum boundary conditions and to perform the simu- 
lations with physical coordinates. In these coordinates 
the cosmological constant, or its generalization as dark 
energy, should be introduced by an explicit term in the 
acceleration equation. The FTM code in its present form 
does not contain a cosmological constant term. This has 
led us to assume the open CDM (OCDM) model with 
Slo = 0.3, /;, = 0.7 and erg = 0.9 (where Hq is the current 
cosmological matter density parameter and erg is the vari- 
ance of the density field convolved with a top- hat window 
of radius 8/i _1 Mpc used to normalize the power spec- 
trum). This is very close to the 'concordance' ACDM 
model in dynamical properties. We are interested here 
in the dynamical evolution of the density profile and its 
dependence on the merging history and therefore the re- 
sults obtained here are valid also in a ACDM cosmology. 
The code was tested in the cosmological context using 
the Santa Barbara Cluster model (Frenk et al. 1999). 

A series of linear constraints on the initial density field 
are used to design the numerical experiments. All the 
constraints are of the same form, namely the value of 
the initial density field at different locations, and evalu- 
ated with different smoothing kernels. We used Gaussian 
kernels for the smoothing procedure, where the width of 
the kernel is fixed so as to encompass a mass M (the 
mass scale on which a constraint is imposed) . The set of 
mass scales and the location at which the constraints are 
imposed define the numerical experiment. 

Assuming a cosmological model and power spectrum of 
the primordial perturbation field, a random realization 
of the field is constructed fro m which a constrained real- 
ization is generated using the IHoffman & Ribakl lfl99lh 
algorithm. Many different realizations of the same exper- 
iment can be performed (Romano-Diaz et al. in prep.). 

3. MODELS 

A set of five different models, i.e., experiments, is de- 
signed here to probe different merging histories of a given 
1O 12 /i _1 M halo in an OCDM cosmology. This halo is 



then constrained to have different substructure on dif- 
ferent mass scales and locations designed to collapse at 
different times. The spherical top-hat model is used here 
to set the numerical value of the constraints. The model 
provides the collapse time of substructures as a function 
of the initial density. This is used only as a general rough 
guide as the various substructures are neither spherical 
nor isolated. Furthermore, the few constraints used here 
do not fully control the experiments. The nonlinear dy- 
namics can in principle affect the evolution in a way not 
fully anticipated from the initial conditions. Even more 
important is the role of the random component of the 
constrained realizations ijHoffman fc Ribak"lll991l) . Thus 
depending on the nature of the constraints and the power 
spectrum assumed, the random component can provide 
other significant substructures at different locations and 
mass scales. This can be handled by adding more con- 
straints and varying their numerical values. 

Our models are designed as follows: Model A (our 
benchmark model) is based on two constraints. One is 
that of a 1Q 12 H~ 1 Mq halo at the origin designed to col- 
lapse at z C oii = 2.1. This halo is embedded in a region 
(2nd constraint) corresponding to a mass of 10 13 /i _1 Mq 
in which the over-density is zero — a region correspond- 
ing to an unperturbed Friedmann model. This is a scale 
larger by about a factor of three (in mass) than the com- 
putational sphere and therefore the constraint cannot be 
exactly fulfilled, yet it constrains the large scale modes 
of the realizations to obey it. These two constraints are 
imposed on all other models. Model B adds two sub- 
structures of mass 5 x 10 11 h~ 1 M & within the 10 12 /i _1 A/ Q 
halo, designed to collapse by z co \\ = 3.7. Model C fur- 
ther splits each ones of the 5 x lO 11 /! -1 !^ halos into 
two 2.5 x lO 11 ^ -1 -/^© substructures (.Zcoii = 5.7). Thus, 
the benchmark halo is design to follow two major merg- 
ers events on its way to virialization. Model D takes the 
Model A and imposes six different small substructures of 
mass 10 11 /i^ 1 M Q scattered within the big halo, designed 
to collapse at about z co \\ = 7.0. Model E attempts to 
simulate a more monolithic collapse in which a nested 
set of constraints, located at the origin, is set on a range 
of mass scales down to M = 10 10 h- 1 M Q (z coll = 6.9). 
All models have been constructed with the same seed of 
the random field. All the density constraints constitute 
(2.5 — 3.5)ct perturbations (where a 2 is the variance of 
the appropriately smoothed field), and were imposed on 
a cubic grid of 128 grid-cells per dimension (Romano- 
Diaz et al., in prep.). Wc evolve the linear initial density 
field from z = 120 until z — by means of the FTM 
code. Since we want to follow as close as possible the 
merging history of our models, we have sampled the sys- 
tem's dynamical evolution with 165 time outputs spaced 
logarithmically in the expansion parameter a. Each halo 
is resolved at z = with around 1.2 x 10 6 particles within 
the virial radius. 

4. RESULTS 

All models differ substantially at early epochs, with 
a subsequent convergence in some of their properties 
but not in others. They lead to the formation of a 
single object of mass ~ 10 12 /i~ 1 M Q via mergers with 
the surrounding substructure and with a slow accretion. 
To analyze the clumpy substructure we define the DM 
halo(s) as having the mean density equal to some crit- 
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ical value A c times the critical density of the universe, 
where A c depends on the redshift and the cosmologi- 
cal model. The top-hat model is used to calculate A c (z) 
and the density is calculated within a virial radius (i? v ; r ). 
The halos are identified initially by the HOP algorithm 
ijEisenstein & Hut Ill998|) and approximated by means of 
a radius i? v ir- Comparison of t he HOP halos with t hose of 
the standard FoF halo finder ijDavis et al. II1985I) shows 
good agreement. Given a group catalog, a merger tree is 
constructed for each model using all the snapshots. 

The evolution toward the final halo is studied by track- 
ing back in time the main branch that leads to this halo. 
The general picture which emerges is that of a halo un- 
dergoing phases of slow and ordered evolution intermit- 
ted by episodes of rapid mass growth via collapse and 
majo r mergers, in agreement with previous sim ulations 
(e.g.. iWechsler et al. ll200l iZhao etalH 12003]) . These 
are referred to as the quiescent and violent phases. 

The structure of a halo is studied by fitting it to an 
NFW profile (Eq. ^) and following the cosmological evo- 
lution of the NFW parameters, namely i? s ,i? v ir- Note, 
that within a given cosmology and per a given halo of 
a given mass, only R s is a free parameter to be fitted. 
The fitting algorithm used here is based on logarithmic 
binning of the halo into spherical shells, and estimating 
R s by minimal % 2 , where the residual in a given shell is 
normalized by its own density. The fitting is performed 
within min(0.6i? V ir, 0.5cZh), where c?h is the distance to 
the nearest massive halo. Although the spherical sym- 
metry of a NFW model ignores some of the dynamical 
properties of a halo, we use it as the first approximation 
to the halo structure. The NFW profile is found to be a 
very good fit to the spherically-symmetric density in the 
quasi-static phases. During the violent phases, such as 
mergers of two almost equal mass halos or collapse of a 
few substructures to form a single halo, the halos are out 
of equilibrium and the NFW fit is only approximate. 

The cosmological evolution of R s and Ry 1T of the 
main halo for all the models is presented in Fig. 1. The 
i? v i r trajectories (i.e., the accretion trajectories) show a 
regular behavior and a linear growth, separated by sud- 
den increases, all of which are associated with the violent 
phases. Most strikingly, R s remains constant in the qui- 
escent phases and grows discontinuously in the violent 
ones. The R s is subject to a ~ 5% — 10% jitter in the 
quiescent phases which grows stronger during the vio- 
lent phases when the halos get out of a dynamical equi- 
librium. The scale density p s shows a similar behavior, 
but in the opposite sense — remaining constant in the 
quasi-static phases and decreasing abruptly in the vio- 
lent ones. We note that all models, except B, converge 
to the same value of R s (within the jitter) at the present 
epoch of a ~ 0.8 — 1, in spite of the different tracks lead- 
ing to it. In Model B, the two major halos have already 
turned around but have yet to merge. Even before this 
last merger, R s is larger than in other models and is ex- 
pected to grow further. 

The present simulations show that there is no clear 
relation between p s and the formation time of the halo, as 
given by the time of the last violent event. This contrasts 
with previous claims about such a correlation. The order 
of models given by the formation time of their halos is 
B, E, A, D, and C. All the models have similar p s except 
for Model C which has a value twice as large (Fig. 1). 



The evolution of the dynamical state of the main halos 
has been analyzed in terms of the internal kinetic en- 
ergy (K) within R s and i? v ir- Any perturbation in the 
halo's internal state (violent mass aggregation, random 
energy acquisition, etc.) will be reflected in its K be- 
havior. We find that K within R s behaves similarly to 
R s — remaining constant during the quasi-static stages 
(within the jitter) and growing discontinuously at the 
violent phases. This implies a possible direct relation 
between the internal kinetic energy and R s . This has 
been tested by comparing the relative changes of K and 
Rs (Figure 2), which shows that the larger is the change 
in K the more R s increases. (Further analysis of this re- 
lation is to be presented elsewhere.) The internal kinetic 
energy computed within i? v i r shows a similar behavior 
but grows very slowly during the quasi-static phases. 

5. DISCUSSION 

The halo growth can be divided into the violent and 
quiescent phases analyzed in our superior time sampled 
A-body simulations using constraintcd realizations. This 
allowed us to show conclusively the details which have 
been hinted about in the literature so far. In this Letter 
we focus on the evolution of R s and i? v ir leaving more 
comprehensive analysis for elsewhere (Romano-Diaz et 
al, in prep.). First, we find that the NFW scale R s 
stays constant during the quiescent phase and changes 
abruptly during the violent one. In contrast, i? v ; r is 
growing linearly in the quiescent and abruptly during the 
violent phases. Second, p s stays unchanged during the 
quiescent phase and drops abruptly during the violent 
phase. Third, the value of R s reflects the violent merg- 
ing history of the halo, and depends on the number of 
violent events and their fractional magnitudes, indepen- 
dent of the time and order of these events. The corollary 
is that p s does not reflect the formation time of the halo. 
Fourth, the fractional change in R s is a nonlinear func- 
tion of the fractional absorbed kinetic energy within R s in 
a violent event. We note, that the accretion trajectories 
in all models converge to the same value. This is a reflec- 
tion of the large-scale structure shared by all the models 
and imposed by the constrained initial conditions. 

The advantage of the constrained realizations lies in 
the unique ability to generate a series of models with one 
or more parameters varied in a controlled manner, while 
all others are held fixed. It allows us to cleanly sepa- 
rate the causc-and-cffect relationship between the initial 
parameters and the outcome of the dynamical evolution. 
This complements the prevailing method of large-scale 
cosmological simulations in which issues of structure and 
evolution arc addressed statistically. The scaling rela- 
tions found in such a statistical analysis are not neces- 
sarily applicable to an individual halo. Here we focus on 
the role of the merging history in the halo evolution, by 
imposing density constraints based on the top-hat model. 
We find that the actual history of a halo does not always 
follow closely the expectations based on the simple top- 
hat model. Inspite of this the different models provide us 
with a good laboratory for experimenting with the role 
of the merging history in shaping the structure of DM 
halos. 

The analogy between the halo evolution and thcrmody- 
namical processes has not escaped our attention. Equat- 
ing the quiescent phases with adiabatic processes and the 
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Fig. 1. — Virial and scale radii behavior (continuous and dotted lines respectively) as function of a for models A, B, C, D and E. The 
discontinuous growths in R a and R v i T match the violent phases that each halo passes through. The horizontal bars represent the mean 
value of R s within the quiescent phases. The square brackets delineate the violent phases. The bottom right panel shows the evolution of 
p B with colors corresponding to the models in other panels. 
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Fig. 2. — The fractional variations in the internal kinetic energy 
SK/K vs those in the scaling radius SR B /R B before and after the 
violent phases. The colors correspond to different generations of 
major mergers. 



violent with non-adiabatic ones leads one to associate 
the behavior of i? s with that of the entropy. In this ter- 
minology the entropy remains constant in the quiescent 
phase and grows discontinuously in the violent phase (see 
Fig. 2). Also, the accretion trajectories play the role of 
adiabats and the system jumps from one adiabat to the 
other by a violent event, not unlike a shock wave. The 
question whether this is just a simple analogy or provides 
a deeper understanding will be addressed elsewhere. 

The results obtained here pertain to the one halo stud- 
ied in the framework of an OCDM cosmology. Yet, 
the conclusions reached from the set of experiments pre- 
sented here are relevant to the understanding of halo for- 
mation in the general CDM cosmologies and in particular 
to the 'benchmark' ACDM cosmology. 
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